1 Introduction

This document serves as a comprehensive guide to replicate the analysis conducted in our experiments. The replication package enables anyone to download, rerun the analysis, and verify our assumptions.

1.1 Data Loading

To begin, load the dataset available in the Data/ directory of our GitHub repository. The dataset is read into R using the following code:

data <- read_sav("Untitled26.sav")

In our dataset: Ranking Variable: Presented by the M201_ prefix, this variable assesses preferences across ten options. Effectiveness Variable: Presented by the M209_ prefix, this variable measures the effectiveness of different options. Easiness Variable: Presented by the M210_ prefix, this variable evaluates the perceived ease of each option.

For more detailed information on these variables and their corresponding options, please refer to our manuscript or questionnaire.

2 Model Design

2.1 Effectiveness

2.1.1 Convert Categorical Variables to Factors

We convert specific columns to ordered factors to ensure correct ordering in categorical analysis.

data$D104 <- factor(data$D104, levels = c("1", "2", "3"), ordered = TRUE)
data$D102 <- factor(data$D102, levels = c("1", "2", "3", "4"), ordered = TRUE)

2.1.2 Data Cleaning - Handle Missing Values

We handle missing values by recoding -9 to NA for the variables starting with “M201_,”M209_“, and”M210_“.

  modeldata <- data %>%
  mutate(across(starts_with("M201_"), ~if_else(. == -9, NA_real_, .))) %>%
  mutate(across(starts_with("M209_"), ~if_else(. == -9, NA_real_, .))) %>%
  mutate(across(starts_with("M210_"), ~if_else(. == -9, NA_real_, .)))

2.1.3 Create New Variables Based on Existing Ones

We create new variables based on existing ones to simplify and categorize the data for analysis.

modeldata <- modeldata %>%
  mutate(
    D104 = D104,
    Expertise = case_when(
      D102 == 1 ~ "Human Factors",
      D102 == 2 ~ "Engineering",
      D102 == 3 ~ "Both",
      TRUE ~ "Other"
    ),
    Dev_Struc = case_when(
      I004_01 >= 1 & I004_01 <= 35 ~ "Waterfall",
      I004_01 >= 36 & I004_01 <= 70 ~ "Hybrid",
      I004_01 >= 71 & I004_01 <= 100 ~ "Agile",
      TRUE ~ "NA"
    ),
    Experience = case_when(
      D104 == 1 ~ "LevelLow",
      D104 == 2 ~ "LevelMedium",
      D104 == 3 ~ "LevelHigh",
      TRUE ~ "Other"
    ),
    AreaOfWork = case_when(
      D101_01 == 2 ~ "Supplier",
      D101_02 == 2 ~ "OEM",
      D101_03 == 2 ~ "Consultancy",
      D101_04 == 2 ~ "Research Institute",
      D101_05 == 2 ~ "University",
      D101_06 == 2 ~ "Government",
      D101_07 == 2 ~ "Non-Automotive",
      TRUE ~ "Other"
    )
  ) %>%
 pivot_longer(cols = starts_with("M209_"), names_to = "placementOption", values_to = "effectivenessRating")

2.1.4 Convert Effectiveness varaible to Factor

We convert the “effectivenessRating” variable to an ordered factor for correct ordering in the analysis.

modeldata$effectivenessRating <- factor(modeldata$effectivenessRating, levels = c("1", "2", "3", "4", "5"), ordered = TRUE)

2.1.5 Define Formula

We define the formula to predict “Effectiveness” based on various predictors.

model_formula <- bf(effectivenessRating ~ placementOption + Experience + Expertise + AreaOfWork + Dev_Struc + (1|ID))

2.1.6 Prior Specification

We specify the prior distributions for the model parameters.

  prior_spec <- c(
  prior(normal(0, 5), "b"),
  prior(normal(0, 10), "Intercept"),
  prior(cauchy(0, 2), "sd")
)

2.1.7 Prior Predictive Checks

We check the and plot the priors from the sample only without empirical data. Sample only from the priors. Evidently the medians are quite evenly set along the \(x\)-axis, and the uncertainty is fairly uniformly distributed among the categories \(1,\ldots,6\) (the bars).

  fit <- brm(
  formula = model_formula,
  data = modeldata,
  family = cumulative(probit), 
  prior = prior_spec,
  chains = 4,
  sample_prior = "only",
  cores = 4,
 seed = 12345
)
## Warning: Rows containing NAs were excluded from the model.
## ld: warning: duplicate -rpath '/Users/amnap/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored
## Start sampling
## Warning: 1 of 4000 (0.0%) transitions ended with a divergence.
## See https://mc-stan.org/misc/warnings for details.
pp_check(fit, type = "bars", nsamples = 250)
## Warning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
## instead.

### Sample with data We fit the model to the data using the specified formula, priors, and settings.

fit <- brm(
  formula = model_formula,
  data = modeldata,
  family = cumulative(probit),
  prior = prior_spec,
  chains = 4,
  cores = 4,
  seed = 12345
)
## Warning: Rows containing NAs were excluded from the model.
## ld: warning: duplicate -rpath '/Users/amnap/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored
## Start sampling

2.1.8 Summary and Diagnostics

We print the summary of the fitted model and perform diagnostic checks to ensure the model has converged and is reliable. Our caterpillar plots look good for all estimated parameters (i.e., they look like fat caterpillars when the chains have mixed well).

mcmc_trace(fit, regex_pars = "^b_") + legend_none()

#Diagnostics such as divergences, tree depth, energy, $\widehat{R}$, and $\mathrm{ESS}$ all look good.

rstan::check_hmc_diagnostics(eval(fit)$fit)
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
if(max(rhat(eval(fit)), na.rm=T) >= 1.01) {
  print("Warning: Rhat >=1.01")
} else {
  print("All Rhat <1.01")
}
## [1] "All Rhat <1.01"
if(min(neff_ratio(eval(fit)), na.rm=T) <= 0.2) {
  print("Warning: ESS <=0.2")
} else {
  print("All ESS >0.2")
}
## [1] "All ESS >0.2"

2.1.9 Posterior Predictive Checks

We perform posterior predictive checks to assess how well the model fits the data.

pp_check(fit, type = "bars", nsamples = 250)
## Warning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
## instead.

2.1.10 Ploting

plot(fit)

## Extract and Plot Posterior Distributions We extract posterior samples and plot the posterior distributions for the model coefficients.

posterior_samples <- posterior_samples(fit)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
coefficients <- grep("b_", colnames(posterior_samples), value = TRUE)
plots <- lapply(coefficients, function(coef_name) {
  ggplot(posterior_samples, aes_string(x = coef_name)) +
    geom_density(fill = "blue", alpha = 0.5) +
    labs(title = paste("Posterior Distribution for", gsub("b_", "", coef_name)),
         x = "Effect Size", y = "Density")
})
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
#for (p in plots) {
 # print(p)
#}

2.1.11 Ploting the posterior probability

Next, let’s plot the posterior probability densities for our population-level estimates on the \(\logit\) scale.

parameter_names <- colnames(posterior_samples)
parameters <- grep("^b_", parameter_names, value = TRUE)
parameters <- parameters[!grepl("Intercept", parameters)]
print(parameters)
##  [1] "b_placementOptionM209_02"      "b_placementOptionM209_03"     
##  [3] "b_placementOptionM209_04"      "b_placementOptionM209_05"     
##  [5] "b_placementOptionM209_06"      "b_placementOptionM209_07"     
##  [7] "b_placementOptionM209_08"      "b_placementOptionM209_09"     
##  [9] "b_placementOptionM209_10"      "b_ExperienceLevelLow"         
## [11] "b_ExperienceLevelMedium"       "b_ExpertiseEngineering"       
## [13] "b_ExpertiseHumanFactors"       "b_ExpertiseOther"             
## [15] "b_AreaOfWorkNonMAutomotive"    "b_AreaOfWorkOEM"              
## [17] "b_AreaOfWorkResearchInstitute" "b_AreaOfWorkSupplier"         
## [19] "b_Dev_StrucHybrid"             "b_Dev_StrucWaterfall"
parameters <- rev(parameters)
p <- mcmc_areas(fit, pars = parameters, prob_outer = 0.95, prob = 0.5) 
#+  vline_0(size = 0.3)
y_labels <- c("Feature Definition", "Feature Realization", "System Design",
              "Dedicated human factors team", "User-oriented teams", "Non-functional requirements team",
              "System/feature evaluation team", "Safety team", "Person/team responsible for the overall system",
              "ExperienceLevelLow", "ExperienceLevelMedium", "ExpertiseEngineering",
              "ExpertiseHumanFactors", "ExpertiseOther", "AreaOfWorkNonMAutomotive",
              "AreaOfWorkOEM", "AreaOfWorkResearchInstitute", "AreaOfWorkSupplier")
y_labels <- rev(y_labels)
p <- p + scale_y_discrete(labels = y_labels) 
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
print(p)

2.1.12 Conditional effects

Below we plot conditional effects for some of our estimates.

ce <- conditional_effects(fit, categorical = TRUE)
p <- conditional_effects(fit, effects = "placementOption", categorical = TRUE)
p_crit <- plot(p)[[1]] +  # Removed plot = FALSE
  scale_color_viridis(discrete = TRUE) +
  scale_fill_viridis(discrete = TRUE) +
  xlab("placementOption") +
  ylab("") +
  scale_x_continuous(breaks=seq(0,1))

ce <- conditional_effects(fit, categorical = TRUE)
p <- conditional_effects(fit, effects = "Experience", categorical = TRUE)
p_crit <- plot(p)[[1]] +  # Removed plot = FALSE
  scale_color_viridis(discrete = TRUE) +
  scale_fill_viridis(discrete = TRUE) +
  xlab("Experience") +
  ylab("") +
  scale_x_continuous(breaks=seq(0,1))

ce <- conditional_effects(fit, categorical = TRUE)
p <- conditional_effects(fit, effects = "Expertise", categorical = TRUE)
p_crit <- plot(p)[[1]] +  # Removed plot = FALSE
  scale_color_viridis(discrete = TRUE) +
  scale_fill_viridis(discrete = TRUE) +
  xlab("Criticality") +
  ylab("") +
  scale_x_continuous(breaks=seq(0,1))

2.2 Easiness

2.2.1 Convert Categorical Variables to Factors

We convert specific columns to ordered factors to ensure correct ordering in categorical analysis.

data$D104 <- factor(data$D104, levels = c("1", "2", "3"), ordered = TRUE)
data$D102 <- factor(data$D102, levels = c("1", "2", "3", "4"), ordered = TRUE)

2.2.2 Data Cleaning - Handle Missing Values

We handle missing values by recoding -9 to NA for the variables starting with “M201_,”M209_“, and”M210_“.

  modeldata <- data %>%
  mutate(across(starts_with("M201_"), ~if_else(. == -9, NA_real_, .))) %>%
  mutate(across(starts_with("M209_"), ~if_else(. == -9, NA_real_, .))) %>%
  mutate(across(starts_with("M210_"), ~if_else(. == -9, NA_real_, .)))

2.2.3 Create New Variables Based on Existing Ones

We create new variables based on existing ones to simplify and categorize the data for analysis.

modeldata <- modeldata %>%
  mutate(
    D104 = D104,
    Expertise = case_when(
      D102 == 1 ~ "Human Factors",
      D102 == 2 ~ "Engineering",
      D102 == 3 ~ "Both",
      TRUE ~ "Other"
    ),
    Dev_Struc = case_when(
      I004_01 >= 1 & I004_01 <= 35 ~ "Waterfall",
      I004_01 >= 36 & I004_01 <= 70 ~ "Hybrid",
      I004_01 >= 71 & I004_01 <= 100 ~ "Agile",
      TRUE ~ "NA"
    ),
    Experience = case_when(
      D104 == 1 ~ "LevelLow",
      D104 == 2 ~ "LevelMedium",
      D104 == 3 ~ "LevelHigh",
      TRUE ~ "Other"
    ),
    AreaOfWork = case_when(
      D101_01 == 2 ~ "Supplier",
      D101_02 == 2 ~ "OEM",
      D101_03 == 2 ~ "Consultancy",
      D101_04 == 2 ~ "Research Institute",
      D101_05 == 2 ~ "University",
      D101_06 == 2 ~ "Government",
      D101_07 == 2 ~ "Non-Automotive",
      TRUE ~ "Other"
    )
  ) %>%
 pivot_longer(cols = starts_with("M210_"), names_to = "placementOption", values_to = "easinessRating")

2.2.4 Convert Easiness varaible to Factor

We convert the “easinessRating” variable to an ordered factor for correct ordering in the analysis.

modeldata$easinessRating <- factor(modeldata$easinessRating, levels = c("1", "2", "3", "4", "5"), ordered = TRUE)

2.2.5 Define Formula

We define the formula to predict “Easiness” based on various predictors.

model_formula <- bf(easinessRating ~ placementOption + Experience + Expertise + AreaOfWork + Dev_Struc + (1|ID))

2.2.6 Prior Specification

We specify the prior distributions for the model parameters.

  prior_spec <- c(
  prior(normal(0, 5), "b"),
  prior(normal(0, 10), "Intercept"),
  prior(cauchy(0, 2), "sd")
)

2.2.7 Prior Predictive Checks

We check the and plot the priors from the sample only without empirical data. Sample only from the priors. Evidently the medians are quite evenly set along the \(x\)-axis, and the uncertainty is fairly uniformly distributed among the categories \(1,\ldots,6\) (the bars).

  fit <- brm(
  formula = model_formula,
  data = modeldata,
  family = cumulative(probit), 
  prior = prior_spec,
  chains = 4,
  sample_prior = "only",
  cores = 4,
 seed = 12345
)
## Warning: Rows containing NAs were excluded from the model.
## ld: warning: duplicate -rpath '/Users/amnap/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored
## Start sampling
## Warning: 1 of 4000 (0.0%) transitions ended with a divergence.
## See https://mc-stan.org/misc/warnings for details.
pp_check(fit, type = "bars", nsamples = 250)
## Warning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
## instead.

### Sample with data We fit the model to the data using the specified formula, priors, and settings.

fit <- brm(
  formula = model_formula,
  data = modeldata,
  family = cumulative(probit),
  prior = prior_spec,
  chains = 4,
  cores = 4,
  seed = 12345
)
## Warning: Rows containing NAs were excluded from the model.
## ld: warning: duplicate -rpath '/Users/amnap/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored
## Start sampling

2.2.8 Summary and Diagnostics

We print the summary of the fitted model and perform diagnostic checks to ensure the model has converged and is reliable. Our caterpillar plots look good for all estimated parameters (i.e., they look like fat caterpillars when the chains have mixed well).

mcmc_trace(fit, regex_pars = "^b_") + legend_none()

#Diagnostics such as divergences, tree depth, energy, $\widehat{R}$, and $\mathrm{ESS}$ all look good.

rstan::check_hmc_diagnostics(eval(fit)$fit)
## 
## Divergences:
## 0 of 4000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 4000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
if(max(rhat(eval(fit)), na.rm=T) >= 1.01) {
  print("Warning: Rhat >=1.01")
} else {
  print("All Rhat <1.01")
}
## [1] "All Rhat <1.01"
if(min(neff_ratio(eval(fit)), na.rm=T) <= 0.2) {
  print("Warning: ESS <=0.2")
} else {
  print("All ESS >0.2")
}
## [1] "Warning: ESS <=0.2"

2.2.9 Posterior Predictive Checks

We perform posterior predictive checks to assess how well the model fits the data.

pp_check(fit, type = "bars", nsamples = 250)
## Warning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
## instead.

2.2.10 Ploting

plot(fit)

## Extract and Plot Posterior Distributions We extract posterior samples and plot the posterior distributions for the model coefficients.

posterior_samples <- posterior_samples(fit)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
coefficients <- grep("b_", colnames(posterior_samples), value = TRUE)
plots <- lapply(coefficients, function(coef_name) {
  ggplot(posterior_samples, aes_string(x = coef_name)) +
    geom_density(fill = "blue", alpha = 0.5) +
    labs(title = paste("Posterior Distribution for", gsub("b_", "", coef_name)),
         x = "Effect Size", y = "Density")
})
#for (p in plots) {
 # print(p)
#}

2.2.11 Ploting the posterior probability

Next, let’s plot the posterior probability densities for our population-level estimates on the \(\logit\) scale.

parameter_names <- colnames(posterior_samples)
parameters <- grep("^b_", parameter_names, value = TRUE)
parameters <- parameters[!grepl("Intercept", parameters)]
print(parameters)
##  [1] "b_placementOptionM210_02"      "b_placementOptionM210_03"     
##  [3] "b_placementOptionM210_04"      "b_placementOptionM210_05"     
##  [5] "b_placementOptionM210_06"      "b_placementOptionM210_07"     
##  [7] "b_placementOptionM210_08"      "b_placementOptionM210_09"     
##  [9] "b_placementOptionM210_10"      "b_ExperienceLevelLow"         
## [11] "b_ExperienceLevelMedium"       "b_ExpertiseEngineering"       
## [13] "b_ExpertiseHumanFactors"       "b_ExpertiseOther"             
## [15] "b_AreaOfWorkNonMAutomotive"    "b_AreaOfWorkOEM"              
## [17] "b_AreaOfWorkResearchInstitute" "b_AreaOfWorkSupplier"         
## [19] "b_Dev_StrucHybrid"             "b_Dev_StrucWaterfall"
parameters <- rev(parameters)
p <- mcmc_areas(fit, pars = parameters, prob_outer = 0.95, prob = 0.5) 
#+  vline_0(size = 0.3)
y_labels <- c("Feature Definition", "Feature Realization", "System Design",
              "Dedicated human factors team", "User-oriented teams", "Non-functional requirements team",
              "System/feature evaluation team", "Safety team", "Person/team responsible for the overall system",
              "ExperienceLevelLow", "ExperienceLevelMedium", "ExpertiseEngineering",
              "ExpertiseHumanFactors", "ExpertiseOther", "AreaOfWorkNonMAutomotive",
              "AreaOfWorkOEM", "AreaOfWorkResearchInstitute", "AreaOfWorkSupplier")
y_labels <- rev(y_labels)
p <- p + scale_y_discrete(labels = y_labels) 
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
print(p)

2.2.12 Conditional effects

Below we plot conditional effects for some of our estimates.

ce <- conditional_effects(fit, categorical = TRUE)
p <- conditional_effects(fit, effects = "placementOption", categorical = TRUE)
p_crit <- plot(p)[[1]] +  # Removed plot = FALSE
  scale_color_viridis(discrete = TRUE) +
  scale_fill_viridis(discrete = TRUE) +
  xlab("placementOption") +
  ylab("") +
  scale_x_continuous(breaks=seq(0,1))

ce <- conditional_effects(fit, categorical = TRUE)
p <- conditional_effects(fit, effects = "Experience", categorical = TRUE)
p_crit <- plot(p)[[1]] +  # Removed plot = FALSE
  scale_color_viridis(discrete = TRUE) +
  scale_fill_viridis(discrete = TRUE) +
  xlab("Experience") +
  ylab("") +
  scale_x_continuous(breaks=seq(0,1))

ce <- conditional_effects(fit, categorical = TRUE)
p <- conditional_effects(fit, effects = "Expertise", categorical = TRUE)
p_crit <- plot(p)[[1]] +  # Removed plot = FALSE
  scale_color_viridis(discrete = TRUE) +
  scale_fill_viridis(discrete = TRUE) +
  xlab("Criticality") +
  ylab("") +
  scale_x_continuous(breaks=seq(0,1))

2.3 Ranking

2.3.1 Convert Categorical Variables to Factors

We convert specific columns to ordered factors to ensure correct ordering in categorical analysis.

data$D104 <- factor(data$D104, levels = c("1", "2", "3"), ordered = TRUE)
data$D102 <- factor(data$D102, levels = c("1", "2", "3", "4"), ordered = TRUE)

2.3.2 Data Cleaning - Handle Missing Values

We handle missing values by recoding -9 to NA for the variables starting with “M201_,”M209_“, and”M210_“.

  modeldata <- data %>%
  mutate(across(starts_with("M201_"), ~if_else(. == -9, NA_real_, .))) %>%
  mutate(across(starts_with("M209_"), ~if_else(. == -9, NA_real_, .))) %>%
  mutate(across(starts_with("M210_"), ~if_else(. == -9, NA_real_, .)))

2.3.3 Create New Variables Based on Existing Ones

We create new variables based on existing ones to simplify and categorize the data for analysis.

modeldata <- modeldata %>%
  mutate(
    D104 = D104,
    Expertise = case_when(
      D102 == 1 ~ "Human Factors",
      D102 == 2 ~ "Engineering",
      D102 == 3 ~ "Both",
      TRUE ~ "Other"
    ),
    Dev_Struc = case_when(
      I004_01 >= 1 & I004_01 <= 35 ~ "Waterfall",
      I004_01 >= 36 & I004_01 <= 70 ~ "Hybrid",
      I004_01 >= 71 & I004_01 <= 100 ~ "Agile",
      TRUE ~ "NA"
    ),
    Experience = case_when(
      D104 == 1 ~ "LevelLow",
      D104 == 2 ~ "LevelMedium",
      D104 == 3 ~ "LevelHigh",
      TRUE ~ "Other"
    ),
    AreaOfWork = case_when(
      D101_01 == 2 ~ "Supplier",
      D101_02 == 2 ~ "OEM",
      D101_03 == 2 ~ "Consultancy",
      D101_04 == 2 ~ "Research Institute",
      D101_05 == 2 ~ "University",
      D101_06 == 2 ~ "Government",
      D101_07 == 2 ~ "Non-Automotive",
      TRUE ~ "Other"
    )
  ) %>%
 pivot_longer(cols = starts_with("M201_"), names_to = "placementOption", values_to = "ranking")

2.3.4 Convert Ranking varaible to Factor

We convert the “ranking” variable to an ordered factor for correct ordering in the analysis.

modeldata$ranking <- factor(modeldata$ranking, levels = c("1", "2", "3", "4", "5", "6","7", "8","9", "10"), ordered = TRUE)

2.3.5 Define Formula

We define the formula to predict “Ranking” based on various predictors.

model_formula <- bf(ranking ~ placementOption + Experience + Expertise + AreaOfWork + Dev_Struc + (1|ID))

2.3.6 Prior Specification

We specify the prior distributions for the model parameters.

  prior_spec <- c(
  prior(normal(0, 5), "b"),
  prior(normal(0, 10), "Intercept"),
  prior(cauchy(0, 2), "sd")
)

2.3.7 Prior Predictive Checks

We check the and plot the priors from the sample only without empirical data. Sample only from the priors. Evidently the medians are quite evenly set along the \(x\)-axis, and the uncertainty is fairly uniformly distributed among the categories \(1,\ldots,6\) (the bars).

  fit <- brm(
  formula = model_formula,
  data = modeldata,
  family = cumulative(probit), 
  prior = prior_spec,
  chains = 4,
  sample_prior = "only",
  cores = 4,
 seed = 12345
)
## Warning: Rows containing NAs were excluded from the model.
## ld: warning: duplicate -rpath '/Users/amnap/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored
## Start sampling
## Warning: 3 of 4000 (0.0%) transitions ended with a divergence.
## See https://mc-stan.org/misc/warnings for details.
pp_check(fit, type = "bars", nsamples = 250)
## Warning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
## instead.

### Sample with data We fit the model to the data using the specified formula, priors, and settings.

fit <- brm(
  formula = model_formula,
  data = modeldata,
  family = cumulative(probit),
  prior = prior_spec,
  chains = 4,
  cores = 4,
  seed = 12345
)
## Warning: Rows containing NAs were excluded from the model.
## ld: warning: duplicate -rpath '/Users/amnap/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored
## Start sampling
## Chain 1 Initialization failed.
## Chain 4 Initialization failed.

2.3.8 Summary and Diagnostics

We print the summary of the fitted model and perform diagnostic checks to ensure the model has converged and is reliable. Our caterpillar plots look good for all estimated parameters (i.e., they look like fat caterpillars when the chains have mixed well).

mcmc_trace(fit, regex_pars = "^b_") + legend_none()

#Diagnostics such as divergences, tree depth, energy, $\widehat{R}$, and $\mathrm{ESS}$ all look good.

rstan::check_hmc_diagnostics(eval(fit)$fit)
## 
## Divergences:
## 0 of 2000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 2000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
if(max(rhat(eval(fit)), na.rm=T) >= 1.01) {
  print("Warning: Rhat >=1.01")
} else {
  print("All Rhat <1.01")
}
## [1] "All Rhat <1.01"
if(min(neff_ratio(eval(fit)), na.rm=T) <= 0.2) {
  print("Warning: ESS <=0.2")
} else {
  print("All ESS >0.2")
}
## [1] "All ESS >0.2"

2.3.9 Posterior Predictive Checks

We perform posterior predictive checks to assess how well the model fits the data.

pp_check(fit, type = "bars", nsamples = 250)
## Warning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
## instead.

2.3.10 Ploting

plot(fit)

## Extract and Plot Posterior Distributions We extract posterior samples and plot the posterior distributions for the model coefficients.

posterior_samples <- posterior_samples(fit)
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
coefficients <- grep("b_", colnames(posterior_samples), value = TRUE)
plots <- lapply(coefficients, function(coef_name) {
  ggplot(posterior_samples, aes_string(x = coef_name)) +
    geom_density(fill = "blue", alpha = 0.5) +
    labs(title = paste("Posterior Distribution for", gsub("b_", "", coef_name)),
         x = "Effect Size", y = "Density")
})
#for (p in plots) {
 # print(p)
#}

2.3.11 Ploting the posterior probability

Next, let’s plot the posterior probability densities for our population-level estimates on the \(\logit\) scale.

parameter_names <- colnames(posterior_samples)
parameters <- grep("^b_", parameter_names, value = TRUE)
parameters <- parameters[!grepl("Intercept", parameters)]
print(parameters)
##  [1] "b_placementOptionM201_02"      "b_placementOptionM201_03"     
##  [3] "b_placementOptionM201_04"      "b_placementOptionM201_05"     
##  [5] "b_placementOptionM201_06"      "b_placementOptionM201_07"     
##  [7] "b_placementOptionM201_08"      "b_placementOptionM201_09"     
##  [9] "b_placementOptionM201_10"      "b_ExperienceLevelLow"         
## [11] "b_ExperienceLevelMedium"       "b_ExpertiseEngineering"       
## [13] "b_ExpertiseHumanFactors"       "b_ExpertiseOther"             
## [15] "b_AreaOfWorkNonMAutomotive"    "b_AreaOfWorkOEM"              
## [17] "b_AreaOfWorkResearchInstitute" "b_AreaOfWorkSupplier"         
## [19] "b_Dev_StrucHybrid"             "b_Dev_StrucWaterfall"
parameters <- rev(parameters)
p <- mcmc_areas(fit, pars = parameters, prob_outer = 0.95, prob = 0.5) 
#+  vline_0(size = 0.3)
y_labels <- c("Feature Definition", "Feature Realization", "System Design",
              "Dedicated human factors team", "User-oriented teams", "Non-functional requirements team",
              "System/feature evaluation team", "Safety team", "Person/team responsible for the overall system",
              "ExperienceLevelLow", "ExperienceLevelMedium", "ExpertiseEngineering",
              "ExpertiseHumanFactors", "ExpertiseOther", "AreaOfWorkNonMAutomotive",
              "AreaOfWorkOEM", "AreaOfWorkResearchInstitute", "AreaOfWorkSupplier")
y_labels <- rev(y_labels)
p <- p + scale_y_discrete(labels = y_labels) 
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
print(p)

2.3.12 Conditional effects

Below we plot conditional effects for some of our estimates.

ce <- conditional_effects(fit, categorical = TRUE)
p <- conditional_effects(fit, effects = "placementOption", categorical = TRUE)
p_crit <- plot(p)[[1]] +  # Removed plot = FALSE
  scale_color_viridis(discrete = TRUE) +
  scale_fill_viridis(discrete = TRUE) +
  xlab("placementOption") +
  ylab("") +
  scale_x_continuous(breaks=seq(0,1))

ce <- conditional_effects(fit, categorical = TRUE)
p <- conditional_effects(fit, effects = "Experience", categorical = TRUE)
p_crit <- plot(p)[[1]] +  # Removed plot = FALSE
  scale_color_viridis(discrete = TRUE) +
  scale_fill_viridis(discrete = TRUE) +
  xlab("Experience") +
  ylab("") +
  scale_x_continuous(breaks=seq(0,1))

ce <- conditional_effects(fit, categorical = TRUE)
p <- conditional_effects(fit, effects = "Expertise", categorical = TRUE)
p_crit <- plot(p)[[1]] +  # Removed plot = FALSE
  scale_color_viridis(discrete = TRUE) +
  scale_fill_viridis(discrete = TRUE) +
  xlab("Criticality") +
  ylab("") +
  scale_x_continuous(breaks=seq(0,1))

3 Model Comparison

Given that our outcome variable follows an ordinal categorical structure, we have numerous modeling options available. Our outcome variable resembles that of an independent achievement model, akin to placing human factors expertise at different positions within an organization. Taking this characteristic into consideration, we anticipate that the following types of models may be suitable options for our analysis (Bürkner and Vuorre, 2019):

  • Cumulative
  • Adjacent-category (acat)
  • Stoping ratio (sratio)
  • Continuation ratio (cratio)
# Fit the model with `cumulative`
fit <- brm(
  formula = model_formula,
  data = modeldata,
  family = cumulative, 
  prior = prior_spec,
  chains = 4,
  cores = 4,
  seed = 12345
)
## Warning: Rows containing NAs were excluded from the model.
## ld: warning: duplicate -rpath '/Users/amnap/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored
## Start sampling
## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 1 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is -93.8941, but should be greater than the previous element, -93.8941 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)
## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 1
## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 1 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is -95.319, but should be greater than the previous element, -95.319 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)
## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 1
## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 1 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is -23.2631, but should be greater than the previous element, -23.2631 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)
## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 1
## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 1 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 2 is -16679.4, but should be greater than the previous element, -16679.4 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)
## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 1
## Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 1 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is -480.442, but should be greater than the previous element, -480.442 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)
## Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 1
## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 2 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is -69.1526, but should be greater than the previous element, -69.1526 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)
## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 2
## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 2 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is -69.421, but should be greater than the previous element, -69.421 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)
## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 2
## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 2 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 4 is -17.6883, but should be greater than the previous element, -17.6883 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)
## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 2
## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 2 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is 6.00503e+258, but should be greater than the previous element, 6.00503e+258 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)
## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 2
## Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 2 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is 7.49579e+06, but should be greater than the previous element, 7.49579e+06 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)
## Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 2
## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 3 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is -81.9327, but should be greater than the previous element, -81.9327 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)
## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 3
## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 3 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is -81.0321, but should be greater than the previous element, -81.0321 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)
## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 3
## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 3 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 5 is -17.7521, but should be greater than the previous element, -17.7521 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)
## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 3
## Chain 3 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 3 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 4 is 6.75391e+17, but should be greater than the previous element, 6.75391e+17 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)
## Chain 3 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 3 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 3
## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 4 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 7 is -28.5643, but should be greater than the previous element, -28.5643 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)
## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 4
## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 4 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 7 is -26.8515, but should be greater than the previous element, -26.8515 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)
## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 4
## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 4 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 3 is inf, but should be greater than the previous element, inf (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)
## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 4
## Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
## Chain 4 Exception: ordered_logistic: Cut-points is not a valid ordered vector. The element at 4 is 5.97504e+26, but should be greater than the previous element, 5.97504e+26 (in '/var/folders/4_/x_f3cbw17fg32f65lnfbvwb0jbn9rd/T/RtmpwZsulx/model-a4b1e71056f.stan', line 100, column 6 to column 63)
## Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
## Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
## Chain 4
print(summary(fit))
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: ranking ~ placementOption + Experience + Expertise + AreaOfWork + Dev_Struc + (1 | ID) 
##    Data: modeldata (Number of observations: 256) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Group-Level Effects: 
## ~ID (Number of levels: 30) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.17      0.14     0.01     0.53 1.00     2080     1952
## 
## Population-Level Effects: 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1]                   -3.00      0.62    -4.21    -1.81 1.00     1785
## Intercept[2]                   -2.12      0.60    -3.29    -0.95 1.00     1707
## Intercept[3]                   -1.48      0.59    -2.63    -0.35 1.00     1682
## Intercept[4]                   -0.95      0.59    -2.09     0.18 1.00     1706
## Intercept[5]                   -0.44      0.58    -1.59     0.67 1.00     1722
## Intercept[6]                    0.07      0.58    -1.07     1.20 1.00     1777
## Intercept[7]                    0.58      0.58    -0.55     1.71 1.00     1786
## Intercept[8]                    1.18      0.58     0.05     2.33 1.00     1800
## Intercept[9]                    2.16      0.60     1.02     3.34 1.00     1873
## placementOptionM201_02         -1.89      0.53    -2.90    -0.87 1.00     1898
## placementOptionM201_03         -0.61      0.52    -1.66     0.37 1.00     1867
## placementOptionM201_04          1.14      0.55     0.04     2.23 1.00     2178
## placementOptionM201_05         -0.98      0.53    -2.03     0.02 1.00     1899
## placementOptionM201_06         -1.47      0.51    -2.48    -0.46 1.00     1796
## placementOptionM201_07          1.05      0.54    -0.00     2.07 1.00     1900
## placementOptionM201_08          0.12      0.50    -0.87     1.11 1.00     1902
## placementOptionM201_09         -0.82      0.51    -1.81     0.17 1.00     1780
## placementOptionM201_10         -1.05      0.52    -2.05    -0.03 1.00     2005
## ExperienceLevelLow             -0.32      0.51    -1.29     0.71 1.00     3582
## ExperienceLevelMedium          -0.15      0.28    -0.69     0.41 1.00     3644
## ExpertiseEngineering           -0.02      0.38    -0.75     0.73 1.00     2992
## ExpertiseHumanFactors           0.12      0.34    -0.56     0.78 1.00     3331
## ExpertiseOther                 -0.05      0.51    -1.07     0.97 1.00     4258
## AreaOfWorkNonMAutomotive       -0.46      0.85    -2.20     1.14 1.00     2963
## AreaOfWorkOEM                  -0.11      0.41    -0.93     0.68 1.00     2175
## AreaOfWorkResearchInstitute     0.13      0.53    -0.91     1.19 1.00     2494
## AreaOfWorkSupplier              0.08      0.43    -0.78     0.93 1.00     2404
## Dev_StrucHybrid                -0.10      0.36    -0.82     0.60 1.00     2895
## Dev_StrucWaterfall              0.08      0.36    -0.62     0.79 1.00     3133
##                             Tail_ESS
## Intercept[1]                    2096
## Intercept[2]                    2009
## Intercept[3]                    2005
## Intercept[4]                    2082
## Intercept[5]                    1905
## Intercept[6]                    2094
## Intercept[7]                    2098
## Intercept[8]                    2163
## Intercept[9]                    2286
## placementOptionM201_02          2540
## placementOptionM201_03          2531
## placementOptionM201_04          2690
## placementOptionM201_05          2461
## placementOptionM201_06          2667
## placementOptionM201_07          2190
## placementOptionM201_08          2548
## placementOptionM201_09          2541
## placementOptionM201_10          2741
## ExperienceLevelLow              2891
## ExperienceLevelMedium           2947
## ExpertiseEngineering            2659
## ExpertiseHumanFactors           3052
## ExpertiseOther                  2885
## AreaOfWorkNonMAutomotive        2656
## AreaOfWorkOEM                   2163
## AreaOfWorkResearchInstitute     2605
## AreaOfWorkSupplier              2654
## Dev_StrucHybrid                 2691
## Dev_StrucWaterfall              3017
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Fit the model with `cratio`
fit1 <- brm(
  formula = model_formula,
  data = modeldata,
  family = cratio, 
  prior = prior_spec,
  chains = 4,
  cores = 4,
  seed = 12345
)
## Warning: Rows containing NAs were excluded from the model.
## ld: warning: duplicate -rpath '/Users/amnap/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored
## Start sampling
print(summary(fit))
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: ranking ~ placementOption + Experience + Expertise + AreaOfWork + Dev_Struc + (1 | ID) 
##    Data: modeldata (Number of observations: 256) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Group-Level Effects: 
## ~ID (Number of levels: 30) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.17      0.14     0.01     0.53 1.00     2080     1952
## 
## Population-Level Effects: 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1]                   -3.00      0.62    -4.21    -1.81 1.00     1785
## Intercept[2]                   -2.12      0.60    -3.29    -0.95 1.00     1707
## Intercept[3]                   -1.48      0.59    -2.63    -0.35 1.00     1682
## Intercept[4]                   -0.95      0.59    -2.09     0.18 1.00     1706
## Intercept[5]                   -0.44      0.58    -1.59     0.67 1.00     1722
## Intercept[6]                    0.07      0.58    -1.07     1.20 1.00     1777
## Intercept[7]                    0.58      0.58    -0.55     1.71 1.00     1786
## Intercept[8]                    1.18      0.58     0.05     2.33 1.00     1800
## Intercept[9]                    2.16      0.60     1.02     3.34 1.00     1873
## placementOptionM201_02         -1.89      0.53    -2.90    -0.87 1.00     1898
## placementOptionM201_03         -0.61      0.52    -1.66     0.37 1.00     1867
## placementOptionM201_04          1.14      0.55     0.04     2.23 1.00     2178
## placementOptionM201_05         -0.98      0.53    -2.03     0.02 1.00     1899
## placementOptionM201_06         -1.47      0.51    -2.48    -0.46 1.00     1796
## placementOptionM201_07          1.05      0.54    -0.00     2.07 1.00     1900
## placementOptionM201_08          0.12      0.50    -0.87     1.11 1.00     1902
## placementOptionM201_09         -0.82      0.51    -1.81     0.17 1.00     1780
## placementOptionM201_10         -1.05      0.52    -2.05    -0.03 1.00     2005
## ExperienceLevelLow             -0.32      0.51    -1.29     0.71 1.00     3582
## ExperienceLevelMedium          -0.15      0.28    -0.69     0.41 1.00     3644
## ExpertiseEngineering           -0.02      0.38    -0.75     0.73 1.00     2992
## ExpertiseHumanFactors           0.12      0.34    -0.56     0.78 1.00     3331
## ExpertiseOther                 -0.05      0.51    -1.07     0.97 1.00     4258
## AreaOfWorkNonMAutomotive       -0.46      0.85    -2.20     1.14 1.00     2963
## AreaOfWorkOEM                  -0.11      0.41    -0.93     0.68 1.00     2175
## AreaOfWorkResearchInstitute     0.13      0.53    -0.91     1.19 1.00     2494
## AreaOfWorkSupplier              0.08      0.43    -0.78     0.93 1.00     2404
## Dev_StrucHybrid                -0.10      0.36    -0.82     0.60 1.00     2895
## Dev_StrucWaterfall              0.08      0.36    -0.62     0.79 1.00     3133
##                             Tail_ESS
## Intercept[1]                    2096
## Intercept[2]                    2009
## Intercept[3]                    2005
## Intercept[4]                    2082
## Intercept[5]                    1905
## Intercept[6]                    2094
## Intercept[7]                    2098
## Intercept[8]                    2163
## Intercept[9]                    2286
## placementOptionM201_02          2540
## placementOptionM201_03          2531
## placementOptionM201_04          2690
## placementOptionM201_05          2461
## placementOptionM201_06          2667
## placementOptionM201_07          2190
## placementOptionM201_08          2548
## placementOptionM201_09          2541
## placementOptionM201_10          2741
## ExperienceLevelLow              2891
## ExperienceLevelMedium           2947
## ExpertiseEngineering            2659
## ExpertiseHumanFactors           3052
## ExpertiseOther                  2885
## AreaOfWorkNonMAutomotive        2656
## AreaOfWorkOEM                   2163
## AreaOfWorkResearchInstitute     2605
## AreaOfWorkSupplier              2654
## Dev_StrucHybrid                 2691
## Dev_StrucWaterfall              3017
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Fit the model with `sratio`
fit2 <- brm(
  formula = model_formula,
  data = modeldata,
  family = sratio, 
  prior = prior_spec,
  chains = 4,
  cores = 4,
  seed = 12345
)
## Warning: Rows containing NAs were excluded from the model.
## ld: warning: duplicate -rpath '/Users/amnap/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored
## 
## Start sampling
print(summary(fit))
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: ranking ~ placementOption + Experience + Expertise + AreaOfWork + Dev_Struc + (1 | ID) 
##    Data: modeldata (Number of observations: 256) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Group-Level Effects: 
## ~ID (Number of levels: 30) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.17      0.14     0.01     0.53 1.00     2080     1952
## 
## Population-Level Effects: 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1]                   -3.00      0.62    -4.21    -1.81 1.00     1785
## Intercept[2]                   -2.12      0.60    -3.29    -0.95 1.00     1707
## Intercept[3]                   -1.48      0.59    -2.63    -0.35 1.00     1682
## Intercept[4]                   -0.95      0.59    -2.09     0.18 1.00     1706
## Intercept[5]                   -0.44      0.58    -1.59     0.67 1.00     1722
## Intercept[6]                    0.07      0.58    -1.07     1.20 1.00     1777
## Intercept[7]                    0.58      0.58    -0.55     1.71 1.00     1786
## Intercept[8]                    1.18      0.58     0.05     2.33 1.00     1800
## Intercept[9]                    2.16      0.60     1.02     3.34 1.00     1873
## placementOptionM201_02         -1.89      0.53    -2.90    -0.87 1.00     1898
## placementOptionM201_03         -0.61      0.52    -1.66     0.37 1.00     1867
## placementOptionM201_04          1.14      0.55     0.04     2.23 1.00     2178
## placementOptionM201_05         -0.98      0.53    -2.03     0.02 1.00     1899
## placementOptionM201_06         -1.47      0.51    -2.48    -0.46 1.00     1796
## placementOptionM201_07          1.05      0.54    -0.00     2.07 1.00     1900
## placementOptionM201_08          0.12      0.50    -0.87     1.11 1.00     1902
## placementOptionM201_09         -0.82      0.51    -1.81     0.17 1.00     1780
## placementOptionM201_10         -1.05      0.52    -2.05    -0.03 1.00     2005
## ExperienceLevelLow             -0.32      0.51    -1.29     0.71 1.00     3582
## ExperienceLevelMedium          -0.15      0.28    -0.69     0.41 1.00     3644
## ExpertiseEngineering           -0.02      0.38    -0.75     0.73 1.00     2992
## ExpertiseHumanFactors           0.12      0.34    -0.56     0.78 1.00     3331
## ExpertiseOther                 -0.05      0.51    -1.07     0.97 1.00     4258
## AreaOfWorkNonMAutomotive       -0.46      0.85    -2.20     1.14 1.00     2963
## AreaOfWorkOEM                  -0.11      0.41    -0.93     0.68 1.00     2175
## AreaOfWorkResearchInstitute     0.13      0.53    -0.91     1.19 1.00     2494
## AreaOfWorkSupplier              0.08      0.43    -0.78     0.93 1.00     2404
## Dev_StrucHybrid                -0.10      0.36    -0.82     0.60 1.00     2895
## Dev_StrucWaterfall              0.08      0.36    -0.62     0.79 1.00     3133
##                             Tail_ESS
## Intercept[1]                    2096
## Intercept[2]                    2009
## Intercept[3]                    2005
## Intercept[4]                    2082
## Intercept[5]                    1905
## Intercept[6]                    2094
## Intercept[7]                    2098
## Intercept[8]                    2163
## Intercept[9]                    2286
## placementOptionM201_02          2540
## placementOptionM201_03          2531
## placementOptionM201_04          2690
## placementOptionM201_05          2461
## placementOptionM201_06          2667
## placementOptionM201_07          2190
## placementOptionM201_08          2548
## placementOptionM201_09          2541
## placementOptionM201_10          2741
## ExperienceLevelLow              2891
## ExperienceLevelMedium           2947
## ExpertiseEngineering            2659
## ExpertiseHumanFactors           3052
## ExpertiseOther                  2885
## AreaOfWorkNonMAutomotive        2656
## AreaOfWorkOEM                   2163
## AreaOfWorkResearchInstitute     2605
## AreaOfWorkSupplier              2654
## Dev_StrucHybrid                 2691
## Dev_StrucWaterfall              3017
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Fit the model wuth `acat`
fit3 <- brm(
  formula = model_formula,
  data = modeldata,
  family = acat, 
  prior = prior_spec,
  chains = 4,
  cores = 4,
  seed = 12345
)
## Warning: Rows containing NAs were excluded from the model.
## ld: warning: duplicate -rpath '/Users/amnap/.cmdstan/cmdstan-2.32.2/stan/lib/stan_math/lib/tbb' ignored
## 
## Start sampling
print(summary(fit))
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: ranking ~ placementOption + Experience + Expertise + AreaOfWork + Dev_Struc + (1 | ID) 
##    Data: modeldata (Number of observations: 256) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Group-Level Effects: 
## ~ID (Number of levels: 30) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.17      0.14     0.01     0.53 1.00     2080     1952
## 
## Population-Level Effects: 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1]                   -3.00      0.62    -4.21    -1.81 1.00     1785
## Intercept[2]                   -2.12      0.60    -3.29    -0.95 1.00     1707
## Intercept[3]                   -1.48      0.59    -2.63    -0.35 1.00     1682
## Intercept[4]                   -0.95      0.59    -2.09     0.18 1.00     1706
## Intercept[5]                   -0.44      0.58    -1.59     0.67 1.00     1722
## Intercept[6]                    0.07      0.58    -1.07     1.20 1.00     1777
## Intercept[7]                    0.58      0.58    -0.55     1.71 1.00     1786
## Intercept[8]                    1.18      0.58     0.05     2.33 1.00     1800
## Intercept[9]                    2.16      0.60     1.02     3.34 1.00     1873
## placementOptionM201_02         -1.89      0.53    -2.90    -0.87 1.00     1898
## placementOptionM201_03         -0.61      0.52    -1.66     0.37 1.00     1867
## placementOptionM201_04          1.14      0.55     0.04     2.23 1.00     2178
## placementOptionM201_05         -0.98      0.53    -2.03     0.02 1.00     1899
## placementOptionM201_06         -1.47      0.51    -2.48    -0.46 1.00     1796
## placementOptionM201_07          1.05      0.54    -0.00     2.07 1.00     1900
## placementOptionM201_08          0.12      0.50    -0.87     1.11 1.00     1902
## placementOptionM201_09         -0.82      0.51    -1.81     0.17 1.00     1780
## placementOptionM201_10         -1.05      0.52    -2.05    -0.03 1.00     2005
## ExperienceLevelLow             -0.32      0.51    -1.29     0.71 1.00     3582
## ExperienceLevelMedium          -0.15      0.28    -0.69     0.41 1.00     3644
## ExpertiseEngineering           -0.02      0.38    -0.75     0.73 1.00     2992
## ExpertiseHumanFactors           0.12      0.34    -0.56     0.78 1.00     3331
## ExpertiseOther                 -0.05      0.51    -1.07     0.97 1.00     4258
## AreaOfWorkNonMAutomotive       -0.46      0.85    -2.20     1.14 1.00     2963
## AreaOfWorkOEM                  -0.11      0.41    -0.93     0.68 1.00     2175
## AreaOfWorkResearchInstitute     0.13      0.53    -0.91     1.19 1.00     2494
## AreaOfWorkSupplier              0.08      0.43    -0.78     0.93 1.00     2404
## Dev_StrucHybrid                -0.10      0.36    -0.82     0.60 1.00     2895
## Dev_StrucWaterfall              0.08      0.36    -0.62     0.79 1.00     3133
##                             Tail_ESS
## Intercept[1]                    2096
## Intercept[2]                    2009
## Intercept[3]                    2005
## Intercept[4]                    2082
## Intercept[5]                    1905
## Intercept[6]                    2094
## Intercept[7]                    2098
## Intercept[8]                    2163
## Intercept[9]                    2286
## placementOptionM201_02          2540
## placementOptionM201_03          2531
## placementOptionM201_04          2690
## placementOptionM201_05          2461
## placementOptionM201_06          2667
## placementOptionM201_07          2190
## placementOptionM201_08          2548
## placementOptionM201_09          2541
## placementOptionM201_10          2741
## ExperienceLevelLow              2891
## ExperienceLevelMedium           2947
## ExpertiseEngineering            2659
## ExpertiseHumanFactors           3052
## ExpertiseOther                  2885
## AreaOfWorkNonMAutomotive        2656
## AreaOfWorkOEM                   2163
## AreaOfWorkResearchInstitute     2605
## AreaOfWorkSupplier              2654
## Dev_StrucHybrid                 2691
## Dev_StrucWaterfall              3017
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Next, we’ll use Leave-One-Out Cross-Validation (LOO) to compare the models:

fit <- add_criterion(fit, criterion = "loo")
fit1 <- add_criterion(fit1, criterion = "loo")
fit2 <- add_criterion(fit2, criterion = "loo")
fit3 <- add_criterion(fit3, criterion = "loo")

loo_compare(fit, fit1, fit2, fit3)
##      elpd_diff se_diff
## fit1  0.0       0.0   
## fit2  0.0       0.0   
## fit  -1.2       3.2   
## fit3 -3.3       3.5

The LOO analysis ranks “fit” as the top-performing model. Based on this ranking, it may be preferable to use “fit” for predictive modeling, as it demonstrates the best performance in terms of predictive accuracy according to the LOO comparison results. Note: We’ve demonstrated this process for “Effectiveness” here, but using the same script, you can perform it for other variables as well. We obtained similar results for all other variables.

4 Computational environment

print(sessionInfo(), locale=FALSE)
## R version 4.3.3 (2024-02-29)
## Platform: x86_64-apple-darwin20 (64-bit)
## Running under: macOS Sonoma 14.4.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] ggthemes_5.1.0     viridis_0.6.5      viridisLite_0.4.2  patchwork_1.2.0   
##  [5] rethinking_2.31    cmdstanr_0.5.3     rstan_2.32.5       StanHeaders_2.32.6
##  [9] mcmcplots_0.4.3    coda_0.19-4.1      mice_3.16.0        lubridate_1.9.3   
## [13] forcats_1.0.0      stringr_1.5.1      purrr_1.0.2        readr_2.1.5       
## [17] tidyr_1.3.1        tibble_3.2.1       tidyverse_2.0.0    brms_2.20.4       
## [21] Rcpp_1.0.12        dplyr_1.1.4        ggplot2_3.5.0      gRain_1.4.1       
## [25] gRbase_2.0.1       bnlearn_4.9.1      haven_2.5.4        bayesplot_1.11.1  
## 
## loaded via a namespace (and not attached):
##   [1] tensorA_0.36.2.1     rstudioapi_0.15.0    jsonlite_1.8.8      
##   [4] shape_1.4.6.1        magrittr_2.0.3       jomo_2.7-6          
##   [7] farver_2.1.1         nloptr_2.0.3         rmarkdown_2.26      
##  [10] vctrs_0.6.5          minqa_1.2.6          denstrip_1.5.4      
##  [13] base64enc_0.1-3      htmltools_0.5.7      distributional_0.4.0
##  [16] curl_5.2.1           broom_1.0.5          mitml_0.4-5         
##  [19] sass_0.4.8           bslib_0.6.1          htmlwidgets_1.6.4   
##  [22] plyr_1.8.9           zoo_1.8-12           cachem_1.0.8        
##  [25] sfsmisc_1.1-17       igraph_2.0.2         mime_0.12           
##  [28] lifecycle_1.0.4      iterators_1.0.14     pkgconfig_2.0.3     
##  [31] colourpicker_1.3.0   Matrix_1.6-5         R6_2.5.1            
##  [34] fastmap_1.1.1        shiny_1.8.0          digest_0.6.34       
##  [37] colorspace_2.1-0     ps_1.7.6             rprojroot_2.0.4     
##  [40] crosstalk_1.2.1      labeling_0.4.3       fansi_1.0.6         
##  [43] timechange_0.3.0     abind_1.4-5          compiler_4.3.3      
##  [46] withr_3.0.0          backports_1.4.1      inline_0.3.19       
##  [49] shinystan_2.6.0      highr_0.10           QuickJSR_1.1.3      
##  [52] pkgbuild_1.4.3       pan_1.9              MASS_7.3-60.0.1     
##  [55] gtools_3.9.5         loo_2.7.0            tools_4.3.3         
##  [58] httpuv_1.6.14        threejs_0.3.3        nnet_7.3-19         
##  [61] glue_1.7.0           nlme_3.1-164         promises_1.2.1      
##  [64] grid_4.3.3           checkmate_2.3.1      reshape2_1.4.4      
##  [67] generics_0.1.3       gtable_0.3.4         tzdb_0.4.0          
##  [70] data.table_1.15.2    hms_1.1.3            utf8_1.2.4          
##  [73] foreach_1.5.2        pillar_1.9.0         markdown_1.12       
##  [76] posterior_1.5.0      later_1.3.2          splines_4.3.3       
##  [79] lattice_0.22-5       survival_3.5-8       tidyselect_1.2.0    
##  [82] miniUI_0.1.1.1       knitr_1.45           gridExtra_2.3       
##  [85] V8_4.4.2             bookdown_0.39        stats4_4.3.3        
##  [88] xfun_0.42            bridgesampling_1.1-2 matrixStats_1.2.0   
##  [91] DT_0.32              stringi_1.8.3        yaml_2.3.8          
##  [94] boot_1.3-29          evaluate_0.23        codetools_0.2-19    
##  [97] cli_3.6.2            RcppParallel_5.1.7   rpart_4.1.23        
## [100] shinythemes_1.2.0    xtable_1.8-4         processx_3.8.3      
## [103] munsell_0.5.0        jquerylib_0.1.4      rstantools_2.4.0    
## [106] ellipsis_0.3.2       dygraphs_1.1.1.6     Brobdingnag_1.2-9   
## [109] lme4_1.1-35.1        glmnet_4.1-8         mvtnorm_1.2-4       
## [112] ggridges_0.5.6       scales_1.3.0         xts_0.13.2          
## [115] crayon_1.5.2         rlang_1.1.3          shinyjs_2.1.0

5 References

Bürkner, P.-C., Vuorre, M., 2019. Ordinal regression models in psychology: A tutorial. Advances in Methods and Practices in Psychological Science 2, 77–101. https://doi.org/10.1177/2515245918823199